Colleen B. Bove1*, Katharine Greene1, Sharla Sugierski2, Nicola G. Kriefall1, Alexa K. Huzar1, Annabel M. Hughes1, Koty Sharp3, Nicole D. Fogarty2, Sarah W. Davies1*

1 Department of Biology, Boston University, Boston, MA USA

2 Department of Biology and Marine Biology, Center for Marine Science, University of North Carolina Wilmington, Wilmington, NC USA

3 Department of Biology, Marine Biology, and Environmental Science, Roger Williams University, Bristol, RI, USA

*Corresponding authors: Colleen B. Bove () and Sarah W. Davies ()


Abstract:

Global change is increasing seawater temperatures and decreasing oceanic pH, driving declines of coral reefs globally. Coral ecosystems are also impacted by local stressors, including microplastics, which are ubiquitous on reefs. While the independent effects of these global and local stressors are well-documented, their interactions remain less explored. Here, we examine the independent and combined effects of global change (ocean warming and acidification) and microplastics exposures on gene expression (GE) and microbial community composition in the endangered coral Acropora cervicornis. Nine genotypes were fragmented and maintained in one of four experimental treatments: 1) ambient conditions (ambient seawater, no microplastics; AMB); 2) microplastics treatment (ambient seawater, microplastics; MP); 3) global change conditions (warm and acidic conditions, no microplastics; OAW); and 4) multistressor treatment (warm and acidic conditions with microplastics; OAW+MP) for 21 days after which corals were sampled for genome-wide GE profiling and ITS and 16S metabarcoding. Overall A. cervicornis GE responses to all treatments were subtle; however, corals in the multistressor treatment exhibited the strongest GE responses, and gene ontology enrichment analyses showcased that genes associated with innate immunity were overrepresented in this treatment. 16S analyses showcased that microbiomes were dominated by the bacterial symbiont Aquarickettsia and were stable, suggesting that A. cervicornis exhibited remarkably low variability in community composition, and future work should focus on functional differences across these bacterial communities as well as the influence of viruses in these responses. Overall, results suggest that local stressors present a unique challenge to endangered coral species under global change.


Citation:

Goes here

Repository DOI hosted on Zenodo [DOI


Microbiome assessment (16S and ITS2) for A. cervicornis microplastics project

This repository contains all code and metadata necessary for the 16S and ITS2 analyses. The code used to assess the gene expression profiles in the manuscript can be found on Sarah Davies’ GitHub repository https://github.com/daviessw/Acer_OAW-Microplastics.


16S Assessment

I am largely following the DADA2 tutorial by Benjamin Callahan that can be found here.

Inspecting the data

Check samples to use

##  [1] "A1"  "A10" "A11" "A12" "A2"  "A3"  "A4"  "A5"  "A6"  "A7"  "A8"  "A9" 
## [13] "B1"  "B10" "B11" "B12" "B2"  "B3"  "B4"  "B5"  "B6"  "B7"  "B8"  "B9" 
## [25] "C1"  "C10" "C11" "C12" "C2"  "C3"  "C4"  "C5"  "C6"  "C7"  "C8"  "C9" 
## [37] "G5"  "G6"



View example quality plots

Viewing the quality of 2 of the forward samples

Viewing the quality of the same 2 of the reverse samples



Filter and trim samples

Perform the filtering

You can see in above figures that the quality of reads drop off towards the end, so we need to filter out these low quality reads

##                  reads.in reads.out
## A1_16S_R1.fastq     46968     44461
## A10_16S_R1.fastq    25343     23991
## A11_16S_R1.fastq    27398     26006
## A12_16S_R1.fastq    45822     44152
## A2_16S_R1.fastq     99390     94836
## A3_16S_R1.fastq     77568     73174

After filtering out the low quality reads, we have maintained about 95.5 of the original reads.



Check the filtering

Viewing the quality of the same 2 of the forward samples post filtering and trimming

Viewing the quality of the same 2 of the reverse samples post filtering and trimming



DADA2 sample processing

Error Rates

The error rates for each possible transition (A→C, A→G, …) are shown. Points are the observed error rates for each consensus quality score. The black line shows the estimated error rates after convergence of the machine-learning algorithm. The red line shows the error rates expected under the nominal definition of the Q-score. We want the estimated error rates (black line) to be a good fit to the observed rates (points), and the error rates to drop with increased quality.



Sample Inference

We are now ready to apply the core sample inference algorithm to the filtered and trimmed sequence data.


Forward sample inference

## Sample 1 - 44461 reads in 9515 unique sequences.
## Sample 2 - 23991 reads in 6094 unique sequences.
## Sample 3 - 26006 reads in 6312 unique sequences.
## Sample 4 - 44152 reads in 7218 unique sequences.
## Sample 5 - 94836 reads in 13509 unique sequences.
## Sample 6 - 73174 reads in 9644 unique sequences.
## Sample 7 - 22069 reads in 4917 unique sequences.
## Sample 8 - 76841 reads in 25268 unique sequences.
## Sample 9 - 41897 reads in 7509 unique sequences.
## Sample 10 - 19892 reads in 5018 unique sequences.
## Sample 11 - 48504 reads in 7342 unique sequences.
## Sample 12 - 32048 reads in 7630 unique sequences.
## Sample 13 - 90475 reads in 7748 unique sequences.
## Sample 14 - 24910 reads in 5725 unique sequences.
## Sample 15 - 34782 reads in 5737 unique sequences.
## Sample 16 - 50069 reads in 5167 unique sequences.
## Sample 17 - 76155 reads in 7895 unique sequences.
## Sample 18 - 34374 reads in 5211 unique sequences.
## Sample 19 - 90460 reads in 8319 unique sequences.
## Sample 20 - 23566 reads in 5815 unique sequences.
## Sample 21 - 53651 reads in 5991 unique sequences.
## Sample 22 - 110660 reads in 10438 unique sequences.
## Sample 23 - 14928 reads in 3601 unique sequences.
## Sample 24 - 30116 reads in 6468 unique sequences.
## Sample 25 - 27927 reads in 5388 unique sequences.
## Sample 26 - 73356 reads in 7647 unique sequences.
## Sample 27 - 23594 reads in 6498 unique sequences.
## Sample 28 - 20343 reads in 4504 unique sequences.
## Sample 29 - 25590 reads in 6768 unique sequences.
## Sample 30 - 10144 reads in 3129 unique sequences.
## Sample 31 - 42054 reads in 7901 unique sequences.
## Sample 32 - 25161 reads in 7367 unique sequences.
## Sample 33 - 73293 reads in 8481 unique sequences.
## Sample 34 - 37175 reads in 5029 unique sequences.
## Sample 35 - 30216 reads in 4299 unique sequences.
## Sample 36 - 43608 reads in 5216 unique sequences.
## Sample 37 - 1795 reads in 421 unique sequences.
## Sample 38 - 162 reads in 63 unique sequences.

Inspecting the returned dada-class object for the first forward sample:

## dada-class: object describing DADA2 denoising results
## 144 sequence variants were inferred from 9515 input unique sequences.
## Key parameters: OMEGA_A = 1e-40, OMEGA_C = 1e-40, BAND_SIZE = 16

The DADA2 algorithm inferred 144 true sequence variants from the 9515 unique sequences in the first sample. There is much more to the dada-class return object than this (see help("dada-class") for some info), including multiple diagnostics about the quality of each denoised sequence variant, but that is beyond the scope of an introductory tutorial.


Reverse sample inference

## Sample 1 - 44461 reads in 9175 unique sequences.
## Sample 2 - 23991 reads in 5924 unique sequences.
## Sample 3 - 26006 reads in 6092 unique sequences.
## Sample 4 - 44152 reads in 7108 unique sequences.
## Sample 5 - 94836 reads in 12281 unique sequences.
## Sample 6 - 73174 reads in 8686 unique sequences.
## Sample 7 - 22069 reads in 4797 unique sequences.
## Sample 8 - 76841 reads in 24793 unique sequences.
## Sample 9 - 41897 reads in 7149 unique sequences.
## Sample 10 - 19892 reads in 4939 unique sequences.
## Sample 11 - 48504 reads in 7069 unique sequences.
## Sample 12 - 32048 reads in 7283 unique sequences.
## Sample 13 - 90475 reads in 7578 unique sequences.
## Sample 14 - 24910 reads in 5526 unique sequences.
## Sample 15 - 34782 reads in 5706 unique sequences.
## Sample 16 - 50069 reads in 5134 unique sequences.
## Sample 17 - 76155 reads in 7407 unique sequences.
## Sample 18 - 34374 reads in 4744 unique sequences.
## Sample 19 - 90460 reads in 8211 unique sequences.
## Sample 20 - 23566 reads in 5706 unique sequences.
## Sample 21 - 53651 reads in 5681 unique sequences.
## Sample 22 - 110660 reads in 10194 unique sequences.
## Sample 23 - 14928 reads in 3754 unique sequences.
## Sample 24 - 30116 reads in 6446 unique sequences.
## Sample 25 - 27927 reads in 5233 unique sequences.
## Sample 26 - 73356 reads in 7528 unique sequences.
## Sample 27 - 23594 reads in 6469 unique sequences.
## Sample 28 - 20343 reads in 4401 unique sequences.
## Sample 29 - 25590 reads in 6315 unique sequences.
## Sample 30 - 10144 reads in 3041 unique sequences.
## Sample 31 - 42054 reads in 7772 unique sequences.
## Sample 32 - 25161 reads in 7026 unique sequences.
## Sample 33 - 73293 reads in 8758 unique sequences.
## Sample 34 - 37175 reads in 5235 unique sequences.
## Sample 35 - 30216 reads in 4296 unique sequences.
## Sample 36 - 43608 reads in 5159 unique sequences.
## Sample 37 - 1795 reads in 465 unique sequences.
## Sample 38 - 162 reads in 69 unique sequences.



Merge paired reads

We now merge the forward and reverse reads together to obtain the full denoised sequences. Merging is performed by aligning the denoised forward reads with the reverse-complement of the corresponding denoised reverse reads, and then constructing the merged “contig” sequences. By default, merged sequences are only output if the forward and reverse reads overlap by at least 12 bases, and are identical to each other in the overlap region (but these conditions can be changed via function arguments).

The mergers object is a list of data.frames from each sample. Each data.frame contains the merged sequence, its abundance, and the indices of the forward and reverse sequence variants that were merged. Paired reads that did not exactly overlap were removed by mergePairs, further reducing spurious output.



Construct sequence table

We can now construct an amplicon sequence variant table (ASV) table, a higher-resolution version of the OTU table produced by traditional methods.


## 
##  221  222  223  227  228  239  246  249  250  251  252  253  254  255  256  257 
##    1    2    4    4    1    1    1    1    1   10  105 3259  140    9    5    2

The sequence table is a matrix with rows corresponding to (and named by) the samples, and columns corresponding to (and named by) the sequence variants. This table contains 3546 ASVs.


After viewing the distribution of read lengths, it looks like we have some that fall outside the expected range (244 - 264) so we will go ahead and remove these non target length sequences.

## 
##  246  249  250  251  252  253  254  255  256  257 
##    1    1    1   10  105 3259  140    9    5    2

This updated table now contains 3533 ASVs.



Remove chimeras

The core dada method corrects substitution and indel errors, but chimeras remain. Fortunately, the accuracy of sequence variants after denoising makes identifying chimeric ASVs simpler than when dealing with fuzzy OTUs. Chimeric sequences are identified if they can be exactly reconstructed by combining a left-segment and a right-segment from two more abundant “parent” sequences.

A total of 40 bimeras were identified from the 3493 input sequences, thus retaining 99.8% of sequences.



Track reads through the pipeline

As a final check of our progress, we can look at the number of reads that made it through each step in the pipeline:

Sample ID Sample Genotype Treatment Input Filtered Denoised Forward Denoised Reverse Merged Nonchimera
A1 15 E - 1 15 MP 46968 44461 44392 44369 44242 44230
A10 6 B - 2 6 OAW+MP 25343 23991 23882 23969 23811 23811
A11 5 B - 1 5 OAW 27398 26006 25982 25986 25976 25976
A12 5a E - 1 5a MP 45822 44152 44099 44109 43973 43973
A2 10a A - 1 10a OAW+MP 99390 94836 94653 94671 94353 94205
A3 5a C - 2 5a AMB 77568 73174 73129 73130 73019 73019
A4 15 E - 2 15 OAW 23219 22069 22048 22045 22026 22026
A5 15 C - 1 15 OAW+MP 82624 76841 76349 76438 73826 72033
A6 10a B - 1 10a MP 44719 41897 41807 41866 41495 41495
A7 15 E - 3 15 AMB 21016 19892 19813 19857 19742 19742
A8 6 D - 1 6 AMB 50868 48504 48363 48414 48090 48044
A9 5a C - 1 5a OAW 33629 32048 31998 31994 31891 31859
B1 11 E - 2 11 OAW+MP 93788 90475 90409 90402 90155 90146
B10 5a A - 1 5a OAW+MP 26111 24910 24846 24859 24708 24683
B11 5 D - 1 5 AMB 36048 34782 34765 34754 34710 34666
B12 12 B - 1 12 OAW+MP 51802 50069 49983 50013 49874 49845
B2 12 E - 1 12 AMB 79689 76155 76112 76109 75821 75770
B3 2 E - 1 2 MP 36469 34374 34349 34353 34289 34289
B4 2 B - 1 2 AMB 93851 90460 90421 90323 90201 90201
B5 6 B - 1 6 MP 24767 23566 23506 23538 23364 23328
B6 2a E - 1 2a OAW 56820 53651 53600 53596 53342 53342
B7 10a C - 1 10a OAW 114866 110660 110347 110599 110242 110229
B8 12 E - 2 12 MP 16088 14928 14822 14864 14767 14690
B9 11 E - 1 11 OAW 31587 30116 30045 30016 29943 29901
C1 5 C - 1 5 MP 29051 27927 27900 27896 27851 27851
C10 2 D - 1 2 OAW+MP 75311 73356 73247 73255 72971 72971
C11 6 E - 1 6 OAW 24571 23594 23583 23578 23571 23571
C12 2a C - 1 2a MP 21209 20343 20305 20329 20214 20185
C2 10a D - 1 10a AMB 26776 25590 25551 25561 25491 25491
C3 11 A - 1 11 MP 10752 10144 10095 10134 10026 10026
C4 5 C - 2 5 OAW+MP 43598 42054 42021 42021 41807 41723
C5 2 A - 1 2 OAW 26542 25161 25122 25096 25002 25002
C6 11 B - 1 11 AMB 76572 73293 73263 73255 72994 72994
C7 2a B - 1 2a OAW+MP 38502 37175 37138 37084 36876 36876
C8 12 A - 1 12 OAW 31488 30216 30071 30042 29692 29692
C9 2a D - 1 2a AMB 44886 43608 43583 43597 43543 43525
G5 NA NA blank_control 1918 1795 1792 1792 1792 1792
G6 NA NA blank_control 169 162 153 158 153 153



Assign taxonomy

It is common at this point, especially in 16S/18S/ITS amplicon sequencing, to assign taxonomy to the sequence variants. The DADA2 package provides a native implementation of the naive Bayesian classifier method for this purpose. The assignTaxonomy function takes as input a set of sequences to be classified and a training set of reference sequences with known taxonomy, and outputs taxonomic assignments with at least minBoot bootstrap confidence.

The dada2 package GitHub maintains the most updated versions of the Silva databases., but I downloaded the databases from the associated Zenodo. The versions in this GitHub repository, used here, were last updated on 26 July 2022.

Let’s inspect the taxonomic assignments:

##      Kingdom    Phylum             Class                 Order              
## [1,] "Bacteria" "Proteobacteria"   "Alphaproteobacteria" "Rickettsiales"    
## [2,] "Bacteria" "Proteobacteria"   "Gammaproteobacteria" "Pseudomonadales"  
## [3,] "Bacteria" "Cyanobacteria"    "Cyanobacteriia"      "Chloroplast"      
## [4,] "Bacteria" "Myxococcota"      "Myxococcia"          "Myxococcales"     
## [5,] "Bacteria" "Proteobacteria"   "Gammaproteobacteria" "Burkholderiales"  
## [6,] "Bacteria" "Campylobacterota" "Campylobacteria"     "Campylobacterales"
##      Family             Genus         Species
## [1,] "Fokiniaceae"      "MD3-55"      NA     
## [2,] "Pseudomonadaceae" "Pseudomonas" NA     
## [3,] NA                 NA            NA     
## [4,] "Myxococcaceae"    "P3OB-42"     NA     
## [5,] "Alcaligenaceae"   "Alcaligenes" NA     
## [6,] NA                 NA            NA

Great! We can now save this and hand it off to Phyloseq for further analyses.



Remove contaminations

Remove mitochondria, chloroplasts, and non-bacteria

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 3277 taxa and 38 samples ]
## sample_data() Sample Data:       [ 38 samples by 5 sample variables ]
## tax_table()   Taxonomy Table:    [ 3277 taxa by 7 taxonomic ranks ]

Removal of mitochondira, chloroplasts, and non-bacteria taxa reduced the total number of taxa from 3493 to 3277 (216 total).



Remove neg control contamination

Just 5 out of the 3277 ASVs were classified as contaminants.


## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 3272 taxa and 36 samples ]
## sample_data() Sample Data:       [ 36 samples by 6 sample variables ]
## tax_table()   Taxonomy Table:    [ 3272 taxa by 7 taxonomic ranks ]



Blast NCBI

I have not rerun this portion since only running the pipeline with the 19 GE samples!

Looks good! We have now cleaned up the sample data to remove contamination from non target organisms and those from the negative controls. The cleaned samples averaged 4.257710^{4} counts (+/- 2.50110^{4} SD).



Rarefy

Rarefy cleaned data

Sample A5 is weirdly high so I am going to remove it from here and through downstream analyses.

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 3272 taxa and 36 samples ]
## sample_data() Sample Data:       [ 36 samples by 6 sample variables ]
## tax_table()   Taxonomy Table:    [ 3272 taxa by 8 taxonomic ranks ]
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 3159 taxa and 36 samples ]
## sample_data() Sample Data:       [ 36 samples by 6 sample variables ]
## tax_table()   Taxonomy Table:    [ 3159 taxa by 8 taxonomic ranks ]



Trim underrepresented OUTs

## [1] "samples with counts below z-score -2.5 :"
## character(0)
## [1] "zscores:"
## named numeric(0)
## [1] "OTUs passing frequency cutoff  1e-04 : 260"
## [1] "OTUs with counts in 0.02 of samples:"
## 
## TRUE 
##  260
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 260 taxa and 36 samples ]
## sample_data() Sample Data:       [ 36 samples by 6 sample variables ]
## tax_table()   Taxonomy Table:    [ 260 taxa by 8 taxonomic ranks ]



Rarefy trimmed data

## named numeric(0)

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 260 taxa and 36 samples ]
## sample_data() Sample Data:       [ 36 samples by 6 sample variables ]
## tax_table()   Taxonomy Table:    [ 260 taxa by 8 taxonomic ranks ]



Microbiome diversity analyses

Notes from phyloseq author Visualize alpha-diversity - Should be done on raw, untrimmed dataset


Relative abundance across samples

Samples appear to be dominated by taxa within the phylum Proteobacteria. Next, I have visualized the diversity of bacteria within Proteobacteria only.

Again, there is a clear dominance across samples by a single order (Rickettsiales), regardless of experimental treatment. I have selected just this order to view the different ASVs that make up this order.



Alpha diversity


Simpson diversity


Shannon diversity

STARTHERE


Species richness


Evenness


All alpha diversity metrics



Beta diversity

## Analysis of Variance Table
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq F value Pr(>F)
## Groups     3 0.12556 0.041855  1.2235 0.3172
## Residuals 32 1.09467 0.034209

## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 1500
## 
## adonis2(formula = seq_all ~ Treatment + Genotype, data = samdf_all, permutations = 1500)
##           Df SumOfSqs      R2      F Pr(>F)
## Treatment  3  0.22644 0.08888 1.0554 0.3664
## Genotype   8  0.60489 0.23742 1.0572 0.3844
## Residual  24  1.71644 0.67370              
## Total     35  2.54777 1.00000



Background microbiome

Because samples are all dominated by Aquarickettsia, we wanted to see if the background diversity varied across treatments once we removed those ASVs.


Relative abundance of non-Aquarickettsia



Alpha diversity


Simpson diversity


Shannon diversity


Species richness


Evenness



Beta Diversity

## Analysis of Variance Table
## 
## Response: Distances
##           Df Sum Sq  Mean Sq F value Pr(>F)
## Groups     3 0.0852 0.028400  1.9262 0.1452
## Residuals 32 0.4718 0.014744

## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 1500
## 
## adonis2(formula = seq_noA_all ~ Treatment + Genotype, data = samdf_noA_all, permutations = 1500)
##           Df SumOfSqs      R2      F Pr(>F)
## Treatment  3   0.5289 0.10597 1.2633 0.1239
## Genotype   8   1.1132 0.22302 0.9971 0.4730
## Residual  24   3.3493 0.67101              
## Total     35   4.9915 1.00000



Core vs. accessory microbiome

Core

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 5 taxa and 36 samples ]
## sample_data() Sample Data:       [ 36 samples by 6 sample variables ]
## tax_table()   Taxonomy Table:    [ 5 taxa by 8 taxonomic ranks ]


Core stats

## Analysis of Variance Table
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq F value Pr(>F)
## Groups     3 0.04572 0.015241  0.4342   0.73
## Residuals 32 1.12330 0.035103

## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 1500
## 
## adonis2(formula = seq_core ~ Treatment + Genotype, data = samdf_core, permutations = 1500)
##           Df SumOfSqs      R2      F Pr(>F)
## Treatment  3   0.5986 0.12424 1.5125 0.1899
## Genotype   8   1.0533 0.21861 0.9980 0.4470
## Residual  24   3.1663 0.65715              
## Total     35   4.8182 1.00000
## Analysis of Variance Table
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq F value Pr(>F)
## Groups     3 0.05829 0.019431   1.142  0.347
## Residuals 32 0.54446 0.017014

## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 1500
## 
## adonis2(formula = seq_core_rel ~ Treatment + Genotype, data = samdf_core, permutations = 1500)
##           Df SumOfSqs      R2      F Pr(>F)
## Treatment  3  0.09704 0.10867 1.2045 0.3438
## Genotype   8  0.15141 0.16957 0.7048 0.7015
## Residual  24  0.64450 0.72176              
## Total     35  0.89295 1.00000

Accessory

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 255 taxa and 36 samples ]
## sample_data() Sample Data:       [ 36 samples by 6 sample variables ]
## tax_table()   Taxonomy Table:    [ 255 taxa by 8 taxonomic ranks ]

Accessory stats

## Analysis of Variance Table
## 
## Response: Distances
##           Df   Sum Sq   Mean Sq F value Pr(>F)
## Groups     3 0.018575 0.0061918  0.7947 0.5059
## Residuals 32 0.249329 0.0077915

## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 1500
## 
## adonis2(formula = seq_acc ~ Treatment + Genotype, data = samdf_acc, permutations = 1500)
##           Df SumOfSqs      R2      F  Pr(>F)  
## Treatment  3   0.8625 0.07497 0.9169 0.68887  
## Genotype   8   3.1176 0.27097 1.2429 0.01399 *
## Residual  24   7.5252 0.65406                 
## Total     35  11.5053 1.00000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Analysis of Variance Table
## 
## Response: Distances
##           Df  Sum Sq   Mean Sq F value Pr(>F)
## Groups     3 0.01390 0.0046333  0.3884 0.7621
## Residuals 32 0.38176 0.0119301

## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 1500
## 
## adonis2(formula = seq_acc_rel ~ Treatment + Genotype, data = samdf_acc, permutations = 1500)
##           Df SumOfSqs      R2      F  Pr(>F)  
## Treatment  3   0.7471 0.07059 0.8415 0.86209  
## Genotype   8   2.7342 0.25833 1.1548 0.07662 .
## Residual  24   7.1027 0.67108                 
## Total     35  10.5841 1.00000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1



Differentially abundant taxa

Performing this analysis with ANCOM and a tutorial for it can be found here. I am mostly following Nicola’s code from her Moorea paper.

No differentially abundant taxa were identified!



ITS2 Assessment

Counts per sample

Raw count statistics of ITS2 samples
mean max min sd
169246 613872 4684 159834
Post-SymPortal processing count statistics of ITS2 samples
mean max min sd
32774 121515 730 31722



Abundance by treatment and genotype



Abundance by treatment only



Manuscript Materials

Methods

ITS2 and 16S metabarcoding

DNA was extracted using an RNAqueous kit (ThermoFisher Scientific) as described above, except samples were not subjected to the DNA removal step. Internal transcribed spacer region 2 (ITS2) PCR amplification was performed using SYM_VAR_5.8S2 and SYM_VAR_REV primers (as described in Hume et al. 2013, 2015) using the following PCR profile: 26 cycles of 95℃ for 40 s, 59℃ for 2 min, 72℃ for 1 min and a final extension of 72℃ for 7 min. A negative control was included and it failed to amplify so it was not sequenced. PCR products were cleaned using the GeneJET PCR Purification kit (ThermoFisher Scientific) according to the manufacturer’s instructions. A second PCR was performed to dual-barcode samples before pooling, which was done based on the visualization of band intensity on a 1% agarose gel. 20 μl of the pool was run on a 2% SYBR Green gel for 100 min at 70 V. The target band was excised and placed in 30 μl of Milli-Q water overnight at 4℃. The supernatant concentration was measured using a DeNovix DS-11+ Spectrophotometer.

The V4 region of the 16S rRNA gene was amplified via PCR using Hyb515f and Hyb806R primers (Parada et al. 2016, Apprill et al. 2015) and the following PCR profile: 30 cycles of 95℃ for 40 s, 63℃ for 2 min, 72℃ for 1 min and a final extension of 7 min. Two negative controls using Milli-Q water were also prepared during the 16S library preparations and these samples were also sequenced. The same procedure as described above for the ITS2 samples was then followed. Concentrations of the ITS2 and 16S pools were used to combine the two pools in a 1:3 ratio respectively. Libraries were sequenced on Illumina MiSeq (paired-end 250 bp) at Tufts Genomics Core Facility.

Demultiplexed reads were pre-processed using bbmap (Bushnell, 2014) to split ITS2 and 16S reads based on primers, while tossing reads that included neither primer. Resulting ITS2 reads were then analyzed by submitting paired fastq.gz files directly to SymPortal, which identifies specific sets of defining intragenomic ITS2 sequence variants (DIVs) to define ITS2 type profiles that are indicative of genetically differentiated Symbiodiniaceae taxa (Hume et al., 2019).

16S primers were removed using cutadapt (Martin, 2011), then DADA2 (Callahan et al., 2016) was used to conduct quality filtering and inference of 3493 sequence variants (ASVs) (see Table S2 to track reads lost through filtering). Taxonomy was assigned at 100% sequence identity using the Silva v. 138.1 database (Quast et al., 2013) and National Center for Biotechnology Information’s nucleotide database using blast+ (Camacho et al., 2009). ASVs matching mitochondria, chloroplasts, or non-bacterial kingdoms were removed (216 total) and 5 ASVs were removed based on negative controls as contaminants (Decontam; (Davis et al., 2018)). Cleaned counts were rarefied to 9409 using vegan (CITE) and trimmed using MCM.OUT (CITE) to remove ASVs present in less than 0.01% of counts, resulting in 260 ASVs across all 36 samples.

Beta diversity of the 16S trimmed counts based on treatment and genotype was assessed using a PCoA on Bray–Curtis dissimilarity (Phyloseq (McMurdie and Holmes 2013)) and a PERMANOVA (vegan; (Oksanen et al., 2020)). We then calculated alpha diversity (Shannon index, Simpson’s index, ASV richness, and evenness) of the rarefied counts using Phyloseq (function estimate_richness(); (McMurdie and Holmes 2013)). Finally, to test for differences in background 16S communities, we removed all ASVs from the dominant genus MD3-55 (Candidatus Aquarickettsia rohweri, referred to as “Aquarickettsia”; (Klinges et al., 2019)) and performed the same beta and alpha diversity assessments. All ITS2 and 16S data, figures, and analyses were completed in (version 3.6.3; R Core Team (2020)) and can be found on GitHub (https://github.com/seabove7/acer_microplastics).



Results

Symbiodiniaceae community composition

Raw ITS2 counts before submission to SymPortal ranged from 4684 to 6.1387210^{5}, with a mean of 1.6924610^{5}. After classification by SymPortal, ITS2 counts ranged from 730 – 1.2151510^{5} per sample across all 36 samples with mean 3.277410^{4} counts. All coral individuals hosted Symbiodinium ‘fitti’ (ITS2 type A3) with 44% of individuals (n=16) hosting small background amounts of Breviolum minutum (ITS2 type B2).


Acropora cervicornis fragments are dominated by order Rickettsiales

Cleaned 16S reads before trimming and rarefaction averaged 42577 ± 25010 (±SD) per sample across all A. cervicornis samples with a minimum of 9409 (Sample 11 A - 1) and a maximum of 107941 (Sample 10a C - 1). Fragments were dominated by taxa from the phylum Proteobacteria (Figure 5A), which was dominated by a single ASV (Genus MD3-55; Ca. Aquarickettsia) within the order Rickettsiales (Figure 5B; see Supplemental Figure XX). Beta diversity analyses did not identify any statistical differences across samples based experimental treatment or genotype (Figure 5C and Table SXX). Similarly, alpha diversity was indistinguishable across treatments (Supplemental Figure SXX and Table SXX). After removing Aquarickettsia from all samples to assess background bacterial communities, there were still no differences in alpha and beta diversity across treatments or genotypes (Supplemental Figure SXX and Table SXX).


Figures

Figure 2

Figure 2. Relative abundance of major ITS2 types across each sample grouped by treatment. Light grey represents Symbiodinium fitti (A3) and dark grey represents Breviolum minutum (B2).



Figure 5

Figure 5: Acropora cervicornis bacterial (16S) relative abundance across fragments and experimental treatments. Diversity at the phylum level (A) depicts dominance of all samples by taxa from Proteobacteria and this (B) phylum was dominated by the order Rickettsiales in most samples. Beta diversity was visualized through (C) multivariate ordination plots (PCoA) of between-sample Bray–Curtis dissimilarity of rarefied ASVs of all taxa and (D) distance to treatment centroids. PCoA ellipses depict 95% confidence intervals and p-values indicate significance of treatment and genotype.



Supplemental figures

Supplemental Figure 1

Supplemental Figure 1: Average seawater conditions across all tanks per treatment throughout the experiment (temperature = °C; dissolved oxygen = ppm; salinity = ppt; pH = total scale). Treatment is depicted by color in all panels: light blue = ambient conditions (AMB), dark blue = ocean acidification and warming (OAW), orange = microplastics (MP), and dark red = OAW and MP (OAW+MP).



Supplemental Figure 5

Supplemental Figure 5: Relative abundance of the four identified ASVs within order Rickettsiales across all samples. Samples were dominated by ASV 37, with background presence of the other 3 ASVs.



Supplemental Figure 6

Supplemental Figure 6: Alpha diversity metrics across all samples in all treatments, depicting (A) Simpson index, (B) Shannon index, (C) species richness (ASVs), and (D) Evenness.



Supplemental Figure 6

Supplemental Figure 5: Alpha (A-D) and beta (E-F) diversity metrics across all samples in all treatments once Aquarickettsia ASVs were removed. Plots depict (A) Simpson index, (B) Shannon index, (C) species richness, (D) Evenness, (E) multivariate ordination plots (PCoA) of between-sample Bray–Curtis dissimilarity of rarefied ASVs and (F) distance to treatment centroids. PCoA ellipses depict 95% confidence intervals and p-values indicate significance of treatment and genotype. Treatment is depicted by color in all panels: blue = ambient conditions (AMB), orange = ocean acidification and warming (OAW), red = microplastics (MP), and dark red = OAW and MP (OAW+MP).



Supplemental tables

Table S1

Table S1: Environmental conditions measured in every tank every day throughout the experiment by treatment and specific to each tank (temperature = °C; dissolved oxygen = ppm ; salinity = ppt; pH = total scale).



Table S4

Table S4: Statistical model output of beta diversity assessments of all ASVs based on (A) treatment and (B) genotype.
  1. ANOVA of all ASV dispersion
Df Sum Sq Mean Sq F value Pr(>F)
Treatment
Groups 3 0.13 0.04 1.22 0.32
Residuals 32 1.09 0.03
Genotype
Groups1 8 0.20 0.02 0.65 0.73
Residuals1 27 1.03 0.04
  1. PERMANOVA of all ASVs
Df SumOfSqs R2 F Pr(>F)
Treatment 3 0.23 0.09 1.06 0.37
Genotype 8 0.60 0.24 1.06 0.38
Residual 24 1.72 0.67
Total 35 2.55 1.00



Table S5

Table S5: Statistical model output of alpha diversity assessments of all ASVs.
Estimate Std..Error t.value
Simpson
(Intercept) 0.67 0.10 6.77
TreatmentMP -0.28 0.13 -2.12
TreatmentOAW -0.23 0.13 -1.76
TreatmentOAW+MP -0.09 0.13 -0.71
Shannon
(Intercept) -0.08 0.28 -0.30
TreatmentMP 0.81 0.36 2.27
TreatmentOAW 0.56 0.36 1.58
TreatmentOAW+MP 0.18 0.36 0.50
Richness
(Intercept) 0.01 0.00 11.96
TreatmentMP 0.00 0.00 0.05
TreatmentOAW 0.00 0.00 -0.91
TreatmentOAW+MP 0.00 0.00 -0.97
Evenness
(Intercept) 0.26 0.08 3.29
TreatmentMP 0.21 0.11 2.01
TreatmentOAW 0.19 0.11 1.77
TreatmentOAW+MP 0.09 0.11 0.81



Table S6

Table S6: Non-Aquarickettsia ASV beta (A-B) and alpha (C) diversity across samples.
  1. ANOVA of non Aquarickettsia dispersion
Df Sum Sq Mean Sq F value Pr(>F)
Treatment
Groups 3 0.09 0.03 1.93 0.15
Residuals 32 0.47 0.01
Genotype
Groups1 8 0.20 0.02 0.65 0.73
Residuals1 27 1.03 0.04
  1. PERMANOVA of non Aquarickettsia
Df SumOfSqs R2 F Pr(>F)
Treatment 3 0.53 0.11 1.26 0.12
Genotype 8 1.11 0.22 1.00 0.47
Residual 24 3.35 0.67
Total 35 4.99 1.00
  1. Alpha diversity metrics of non Aquarickettsia
Estimate Std..Error t.value Pvalue
Simpson
(Intercept) 0.12 0.02 7.38 0.00
TreatmentMP 0.01 0.02 0.32 0.75
TreatmentOAW -0.02 0.02 -0.92 0.36
TreatmentOAW+MP -0.01 0.02 -0.76 0.45
Shannon
(Intercept) 1.57 0.03 45.48 0.00
TreatmentMP 0.00 0.05 0.07 0.95
TreatmentOAW 0.05 0.05 1.15 0.26
TreatmentOAW+MP 0.05 0.05 1.06 0.30
Richness
(Intercept) 3.09 0.10 31.45 0.00
TreatmentMP 0.07 0.12 0.61 0.55
TreatmentOAW 0.15 0.12 1.27 0.22
TreatmentOAW+MP 0.15 0.12 1.28 0.21
Evenness
(Intercept) 1.26 0.04 35.25 0.00
TreatmentMP 0.03 0.04 0.56 0.58
TreatmentOAW -0.03 0.04 -0.59 0.55
TreatmentOAW+MP -0.01 0.04 -0.34 0.74



Session Information

All code was written by Colleen B. Bove, feel free to contact with questions.

Session information from the last run date on 24 November 2022:

## R version 3.6.3 (2020-02-29)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS  10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats4    parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] lmerTest_3.1-3              Rmisc_1.5.1                
##  [3] plyr_1.8.7                  nlme_3.1-155               
##  [5] webshot_0.5.2               shades_1.4.0               
##  [7] ggh4x_0.2.2.9000            compositions_2.0-4         
##  [9] patchwork_1.1.1             lme4_1.1-28                
## [11] car_3.1-0                   carData_3.0-5              
## [13] performance_0.8.0           awtools_0.2.1              
## [15] microbiome_1.8.0            ggpubr_0.4.0               
## [17] MCMC.OTU_1.0.10             MCMCglmm_2.33              
## [19] ape_5.6-1                   coda_0.19-4                
## [21] Matrix_1.3-4                janitor_2.1.0              
## [23] vegan_2.5-7                 lattice_0.20-45            
## [25] permute_0.9-7               plotly_4.10.0              
## [27] RColorBrewer_1.1-3          decontam_1.6.0             
## [29] phyloseq_1.30.0             kableExtra_1.3.4.9000      
## [31] ShortRead_1.44.3            GenomicAlignments_1.22.1   
## [33] SummarizedExperiment_1.16.1 DelayedArray_0.12.3        
## [35] matrixStats_0.62.0          Biobase_2.46.0             
## [37] Rsamtools_2.2.3             GenomicRanges_1.38.0       
## [39] GenomeInfoDb_1.22.1         Biostrings_2.54.0          
## [41] XVector_0.26.0              IRanges_2.20.2             
## [43] S4Vectors_0.24.4            BiocParallel_1.20.1        
## [45] BiocGenerics_0.32.0         forcats_0.5.1              
## [47] stringr_1.4.0               dplyr_1.0.8                
## [49] purrr_0.3.4                 readr_2.1.2                
## [51] tidyr_1.2.0                 tibble_3.1.7               
## [53] ggplot2_3.3.6               tidyverse_1.3.1            
## [55] dada2_1.20.0                Rcpp_1.0.9                 
## [57] knitr_1.40                 
## 
## loaded via a namespace (and not attached):
##   [1] utf8_1.2.2             tidyselect_1.1.1       htmlwidgets_1.5.4     
##   [4] grid_3.6.3             Rtsne_0.15             munsell_0.5.0         
##   [7] codetools_0.2-18       ragg_1.1.3             interp_1.0-33         
##  [10] withr_2.5.0            colorspace_2.0-3       highr_0.9             
##  [13] rstudioapi_0.13        robustbase_0.93-9      bayesm_3.1-4          
##  [16] ggsignif_0.6.3         labeling_0.4.2         GenomeInfoDbData_1.2.2
##  [19] hwriter_1.3.2.1        farver_2.1.1           rhdf5_2.30.1          
##  [22] vctrs_0.4.1            generics_0.1.3         xfun_0.29             
##  [25] R6_2.5.1               bitops_1.0-7           cachem_1.0.6          
##  [28] assertthat_0.2.1       scales_1.2.0           gtable_0.3.0          
##  [31] processx_3.5.2         rlang_1.0.4            systemfonts_1.0.2     
##  [34] splines_3.6.3          rstatix_0.7.0          lazyeval_0.2.2        
##  [37] broom_1.0.0            yaml_2.3.5             reshape2_1.4.4        
##  [40] abind_1.4-5            modelr_0.1.8           crosstalk_1.2.0       
##  [43] backports_1.4.1        tensorA_0.36.2         tools_3.6.3           
##  [46] cubature_2.0.4.2       ellipsis_0.3.2         jquerylib_0.1.4       
##  [49] biomformat_1.14.0      zlibbioc_1.32.0        RCurl_1.98-1.7        
##  [52] ps_1.6.0               deldir_1.0-6           cowplot_1.1.1         
##  [55] haven_2.4.3            cluster_2.1.2          fs_1.5.2              
##  [58] magrittr_2.0.3         data.table_1.14.2      reprex_2.0.1          
##  [61] hms_1.1.1              evaluate_0.15          jpeg_0.1-9            
##  [64] readxl_1.3.1           gridExtra_2.3          compiler_3.6.3        
##  [67] crayon_1.5.1           minqa_1.2.4            htmltools_0.5.2       
##  [70] mgcv_1.8-38            corpcor_1.6.10         tzdb_0.2.0            
##  [73] RcppParallel_5.1.5     lubridate_1.8.0        DBI_1.1.3             
##  [76] dbplyr_2.1.1           MASS_7.3-55            boot_1.3-28           
##  [79] ade4_1.7-18            cli_3.3.0              insight_0.18.0        
##  [82] igraph_1.2.11          pkgconfig_2.0.3        numDeriv_2016.8-1.1   
##  [85] xml2_1.3.3             foreach_1.5.2          svglite_2.1.0         
##  [88] bslib_0.4.0            multtest_2.42.0        rvest_1.0.3           
##  [91] snakecase_0.11.0       callr_3.7.0            digest_0.6.29         
##  [94] rmarkdown_2.13         cellranger_1.1.0       nloptr_1.2.2.3        
##  [97] lifecycle_1.0.1        jsonlite_1.7.2         Rhdf5lib_1.8.0        
## [100] viridisLite_0.4.0      fansi_1.0.3            pillar_1.8.0          
## [103] fastmap_1.1.0          httr_1.4.3             DEoptimR_1.0-11       
## [106] survival_3.2-13        glue_1.6.2             png_0.1-7             
## [109] iterators_1.0.14       stringi_1.7.8          sass_0.4.0            
## [112] textshaping_0.3.6      latticeExtra_0.6-30